Metastability and anomalous fixation in evolutionary games on scale-free networks 
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We study the influence of complex spatial structure on the metastability and fixation properties 
of a set of evolutionary processes characterized by frequency-dependent selection. In the frame- 
work of evolutionary game theory, we analyze the dynamics of snowdrift games (characterized by a 
metastable coexistence state) on scale-free networks. Using an effective diffusion theory we demon- 
strate how the complex structure of the network affects the system's metastable state and leads to 
anomalous fixation. In particular, we analytically and numerically show that the probability and 
mean time of fixation are characterized by stretched exponential behaviors with exponents depend- 
ing nontrivially on the network's degree distribution. Our approach is also shown to be applicable 
to models, like coordination games, characterized by the absence of metastability prior to fixation. 
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The dynamics of systems where successful traits spread 
at the expense of others is naturally modeled in the 
framework of evolutionary game theory (EGT) Q. EGT 
involves frequency-dependent selection, where selection 
depends on the species instantaneous concentration. Tra- 
ditionally EGT is studied in terms of differential equa- 
tions, and is an approach transcending almost every as- 
pect of evolutionary biology fl','3|. While the EGT classic 
setting was originally proposed to describe the evolution 
of infinitely large and spatially homogeneous populations, 
it is known that evolutionary dynamics is affected by de- 
mographic noise and by the population's spatial arrange- 
ment (structure) 3-5]. Evolutionary dynamics is often 
characterized by the central notion of fixation, first in- 
troduced in population genetics Q, that refers to the 
possibility that a "mutant type" takes over the entire 
population. In contrast to spatially-homogeneous (well- 
mixed) populations, accounting for the population's spa- 
tial arrangement can give rise to various scenarios: e.g., 
in games modeling social dilemmas, the local interactions 
on regular lattices may enhance or inhibit the resistance 
of "cooperators" against the invasion by "defectors" [Jl . 
In this context, evolutionary dynamics on networks (7l| 
provides a general and unifying framework to describe 
the dynamics of both well- mixed and spatially-structured 
populations [sj. In spite of their importance, fixation 
of evolutionary processes on networks have been mostly 
studied in idealized situations, e.g. for two-state systems 
under a constant weak selective bias In these 

works, it has been shown that the update rules and the 
network structure effectively renormalize the population 
size and affect the fixation properties. While these results 
are of great interest, they do not account for frequency- 
dependent selection, which is an important evolutionary 
mechanism and may lead to a long-lived metastable 
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(or coexistence) state prior to fixation [12|, |13| . As a firm 
understanding of metastability on complex networks is 
still lacking, we here study metastability and fixation on 
a class of scale-free networks in the EGT framework. Our 
findings will also be directly relevant to the dynamics of 
epidemic outbreaks [14] and population genetics 15 1. 

Here, we study "snowdrift games" (SGs, see be- 
low) fll, K3\ that are the paradigmatic EGT models ex- 
hibiting metastability (see [Tq] for their experimental rel- 
evance). In the case of well-mixed populations (complete 
graphs) , the fixation properties of SGs typically exhibit 
an exponential dependence on the population size, see 
e.g. [13|. Using an individual- based approach, we show 
that the spatial structure of scale-free networks yields 
anomalous fixation and metastability characterized by a 
stretched exponential dependence on the population size, 
in stark contrast with their non-spatial counterparts. We 
also show that such a dependence characterizes fixation 
in models like coordination games which do not possess 
a long-lived metastable state prior to fixation [ij . 

The model. We consider a network comprising N 
nodes, each of which is either occupied by an individ- 
ual of type C (cooperator) or by a D-individual (defec- 
tor). The occupancy of the node i is encoded by the 
random variable rji, with rji = 1 if the node i is occupied 
by a C and rji = otherwise. The state of the system 
is thus described by {ry} = {rii}^ and the density of co- 
operators present in the system is p = X^t^i Vi/^- The 
network is specified by its adjacency matrix A = [Aij], 
whose elements are 1 if the nodes ij are connected and 
otherwise. The network is also characterized by its de- 
gree distribution = N^/N, where Nf^ is the number of 
nodes of degree k. EGT is traditionally concerned with 
large and homogeneous populations (i.e. TV — >■ oo and 
Aij ^ l,yij) whose mean field dynamics is described by 
the celebrated rephcator equation [l|, {d/dt)p{t) = 
p{t){l - p{t))[U^{p) - n^(p)], where 11^ /^{p) are the 
cooperator/defector average payoffs derived from the 
game's payoff matrix. For a generic two-strategy cooper- 
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ation dilemma, the payoff of C against another C is de- 
noted a and that of D playing against D is d. When C and 
D play against each other, the former gets payoff h and 
the latter gets c [H- Here, we are particularly interested 
in SGs, for which c > a and b > d. SGs are characterized 
by a stable interior fixed point p.^, = {d~b)/{a — b — c + d) 
and unstable absorbing states p — (all-D) and p = 1 
(all-C). When the population size is finite {N < oo), 
the role of fluctuations is important and becomes a 
metastable state whose decay time on coniplete graphs 
(Aij = l,Vij) grows exponentially with N [13j . 

In a spatial setting, the interactions are among nearest- 
neighbor individuals and the species payoffs are defined 
locally: C and D players at node i interacting with a 
neighbor at site j respectively receive payoffs 11^ = 
arjj + 5(1 — r]j) and = crjj + d{l — rjj). In the spirit of 
the Moran model (in the weak selection limit) 0, 
each species local reproductive potential, or fitness, is 

C / D 

given by the difference of 11^^ relative to the popula- 
tion mean payoff Ily(i). Here, we make the mean-field- 
like choice IL,j{t) = p(i)Hg + {1 - p(i))H,^ to include 
what arguably is the simplest mechanism ensuring the 
formation of metastability. It is customary to introduce 
a selection strength s in the definition of the fitness to 
unravel the interplay between random fluctuations and 
selection 0, d, Here, the fitnesses of C/D at node i 
interacting with a neighbor at site j are 

ffj = 1 + s[Hg - H,,] and = 1 + s[H,^ - (1) 

These expressions comprise a baseline contribution (set 
to 1) and a selection term proportional to the relative 
payoffs. Moreover, we consider a system evolving accord- 
ing to the so-called "link dynamics" (LD) 0, [l3|: a link 
is randomly selected at each time step and if it connects 
a CD pair, one of the neighbors is randomly selected for 
reproduction with a rate proportional to its fitness while 
the other is replaced by the offspring [ITJ . Fixation under 
LD is unaffected by the network topology when the selec- 
tion is frequency- independent [lo| . However, in the more 
generic case of frequency-dependent selection considered 
here, metastability and fixation properties are drastically 
altered by the network's structure, as shown below. 

The evolutionary dynamics of the system will be ana- 
lyzed in terms of the (subgraph) density pk — Yl'iVil^k 
of cooperators C on nodes of degree k (the prime means 
that the sum is restricted to nodes of degree k). Quanti- 
ties necessary for our analysis are the m*''— moment of the 
degree distribution, p.^ = Ylk ^™'^fc = Yli l^i where 
ki denotes the degree of node i, and the degree- weighted 
density of cooperators w = ^^{k/ pi)nkPk- 

Effective diffusion theory. To implement the evolu- 
tionary dynamics, we introduce — rji[\ — rjj)f'j< and 



(1 — 'ili)Vjfiji where 77.^(1 — rjj) is non-zero only 



when the nodes ij are occupied by a CD pair. In the LD, 
the probability to select the neighbor j of node i for an 
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FIG. 1: (Color online). Densities p (thin solid), uj (dashed), 
pi (thick solid), p2 (dotted), ps (dashed-dotted), on a scale- 
free network with = 3 for a SG with a = d = l,b = 1.4, c = 
1.6. Initially Pfc>^i(0) = l,pfc<^i(0) = 0. Top panel: Single 
realization for s — 0.001 and A'' — 10*. Left/right lower 
panels: Trajectories averaged over 500 samples (A^ = 2000) 
for s = 0.001 (top), s = 0.01 (middle), s = 0.1 (bottom). 
Right: p and pi versus oj (the straight line is an eye guide). 



update is Aij / (N pi) and the transition r]i — > l — rji hence 
occurs with probability Y.j 4jk + ^i^l tlfll- The 
subgraph density pk changes hy±6pk — ±1/Nk accord- 
ing to a birth-death process [18| defined by the transi- 
tion rates T+{pk) = J2[ A,,^,,/{Npi) and T-{pk) = 
respectively. For our analytical 
treatment, we focused on degree-heterogeneous networks 
with degree-uncorrelated nodes, as in Molloy-Reed net- 
works (MRN) [Til, yielding A,j = hkj/{Npi). Our nu- 
merical simulations were performed using the "redirec- 
tion algorithm" that generates degree-correlated scale- 
free networks [2^. Yet, it has been shown that the dy- 
namics of the latter is close to that on MRN 10]. With 
N^^ — UkPk, the transition rates become 



T+{pk) = T+ = {uk/pi) [1 + s{b - d){l - p)] k{l - pk)oj 
T^{pk) = JV = {rik/pi) [1 - s(a - c)p] kpk{l - w). (2) 

We notice that are nonzero provided that the mean 
degree pi does not diverge with N — >• 00 [2l|. In the 



limit of weak selection intensity (s ^ 1) ^l5|, one can use 
the diffusion theory to treat the birth-death process de- 
fined by ([2]) . This yields a multivariate backward Fokker- 
Planck equation (FPE) whose generator reads 



in-Tk) d , (T++r,7) d 
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dpk 



2Nnl 



dpi 



, (3) 



with time increments 5t — [l^. Moreover, 

when the selection intensity is weak (s ^ 1), a 
timescale separation allows to greatly simplify the 
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analysis [10|, [llj. Indeed, when t ^ s~ the se- 
lection pressure is negligible and the quantity p is 
conserved [§]. In fact, using ^ at mean field level gives 
{d/dt)p = s(a — 6 — c + (i)w(l — ui)[p — p*) (where the 
upper bar denotes the ensemble average). This indicates 
that p relaxes to its metastable value on a timescale 
see Fig. ^ At mean field level, Eqs. (|2|l 
also yields {d/dt)pk = {T^{pk) -T^{Pk))/nk = (fc/Mi) 
x{Q - pk + s[{b-d)u}[l-p){l~pk) + (a-c)(l-w)pA;p]}. 
This indicates that after a timescale of order 0(1), 
Pk ~ w, and also p « w since p = J2kPk''^>'- With 
Pk ~ u) ~ p, the rate equation for pfe becomes 
{d/dt)pk ~ -(fc/^i)(&-d)s(l-pfc)pfc(pfc/p,-l). Hence, 
while after a time of order 0(1), Pk ~ 'I' ~ p, all these 
quantities approach p^ on a longer timescale t ~ s~^. 
This behavior is corroborated by numerical simulations, 
see Fig. [T] [22]. Because fixation occurs on much longer 
time-scales than (see below), we henceforth use 

the approximation that on average pk ~ p ~ uj. Using 
Eq. ([U, changing variables pk ^ and using the 
definition of lu that yields dp,. — > (knk / fJ-i)duj , we arrive 
at the effective single- coordinate FPE generator 

L0{1 — Lu) 



iVeff 



(4) 



Here cr = 2(6— d)NcsScs/ p*, and the effective popu- 
lation size and selection intensity are given by N^s = 
N {pi)^/p3 and SeS = s p2/{pi)'^- For scale-free net- 
works with degree distribution Uk oc k~'^ and finite av- 
erage degree (i.e. 7 > 2) [iH, the maximum degree is 
kmax ~ iVi/^'^-i) This can be used to calculate the 
moments pm [13 and obtain the scaling of A^cffScg: 



^b-d) ,jPiP2 
a = sN 



p* 



M3 



sN, v> 4 

sN/\nN, v = 4: 

^^(2.-5) /(u-l)^ 3<!.<4(5) 

sVNlnN, i/ = 3 

g^(>.-2)/(i.-i)^ 2 < < 3. 



We have checked analytically and by numerical simula- 
tions that our theory is valid when s'^g <C N~g. Hence, 
for 2 < i/ < 4 it is applicable over a broader range of val- 
ues of s than on complete graphs; e.g. ^ jsf-i/i'^-i) 
when 2 < < 3 (s^ <C on complete graphs [isf). 

Fixation properties. Evolutionary dynamics is char- 
acterized by the fixation probability 4>'~'' (uj) that a sys- 
tem with initial degree-weighted density uj is taken over 
by cooperators. In the framework of the effective diffu- 
sion theory and using (H)) the fixation probability obeys 
5eff(w)(/)*" (w) = with the boundary conditions <j)'"{0) = 
[isj . The solution of this equation is 



0^(w) 



erfi [p^^/a]- erfi [{p^ - uj)^/a] 
erfi [p*y^] + erfi [(1 - p*)y^] ' 



(6) 
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FIG. 2: (Color online). Fixation probability (f versus 
rescaled arc oc a [see ((5])] for SGs with a = d = 1, 6 = 1.05, 
c = 1.075 (p. = 0.4) and s = 0.075 (□), 0.1 (A), 0.125 (o) 
and 0.15 (o). Numerical results for u = 2.5 (top), f — 2.75 
(middle), u = 3 (bottom) collapse linearly (dashed) with 
—sN^ip2/p3, in agreement with the theory, see text. Here 
N = 1000 - 4000 and initially = p = u: = 100/iV. 



where erfi(z) = e" du. Let us consider the (bio- 

logically relevant) case of a small initial densityof co- 
operators such that w ^ 1, weak selection and 
a large population such that pjcr » 1 and metastabil- 
ity is guaranteed. Using the asymptote erfi(a;) ~ for 
X ^ 1 in Eq. ([6]), we distinguish two cases: (i) when 
p* < 1/2, In^*^ ~ -(1 - 2p*)cr; (u) when p* > 1/2 
and Lo > 2p^, — 1, ln(l — (j)'^) ~ —{2p^, — 1)0-, while 
ln(l — (f)'^') ~ — w(2p* — uj)a if p* > 1/2 and oj < 2p* — 1. 
The stretched-exponential dependence of (/j*^ on N is 
shown in Fig. [2] where (as in Figs. [3] and data have 
been rescaled and collapsed with a oc sN pip2/ P3 along 
a single line. Since In^*^ ^ — on complete graphs [l^, 
our results demonstrate how the complex structure of the 
network drastically affects the fixation probability. 

Another quantity of great interest is the (uncondi- 
tional) mean fixation time (MET) r(w) - the mean time 
necessary to reach an absorbing boundary. Here, using 
Eq. dH) the MET is obtained by solving t/eff (w)r(w]_ = -1 



with the boundary conditions t(0) = t(1) = 
Using standard methods 0, [l3| , we obtain 



r(cj) = 2Afeff 



-e(y) ry 

(l-(/)^'(w)) / dy^^ -I dze^'-'^ 



v{^-y)Jo 



+ (j)^{uj) / dy 



3-e(y) 



dz e 



e{z) 



; 6(z) = (Tz(z — 2p*). 



For pjcr 3> 1 the inner integrals can be computed by 
expanding Q{z) around its extremal values (z = for z e 
[0,Ph.] and z = 1 for z S [p*, 1]), while the outer integral 
is computed via the saddle-point approximation around 
Lj = p^. To leading order, one thus obtains t(u)) ~ (1 — 
(f)^ {uj))e'^p'' when u > p^ and t{uj) - 4>'^ {uj)e'^^'^-P''>^ 
otherwise. When a ^ 1 and p* < 1/2 (as in Fig. [3]), we 
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FIG. 4: (Color online). Fixation probability 4> versus arc oc 
(7 for CG with V = 2.5, a = 1.6, h = c = 1 and d = 1.2 
(p, = 0.75) and s = 0.05 (□), 0.075 (A), 0.1 (o) and 0.125 (o). 
Numerics are consistent with the theoretical result In^*^ ~ 
-siV^i/.i2/At3. Here = 1000 - 4000 and pfe(O) = 0.5. 



FIG. 3: (Color online). MFT r versus rescaled are oc a 
for V = 2.5 (top), v = 2.75 (middle), and = 3 (bottom). 
Numerical results (symbols) collapse along lines (dashed) in 
agreement with 0, see text. Here, symbols and parameters 
are the same as in Fig. [2l and initially pk = p = uj = 0.5. 



structure on the evolutionary dynamics. 



find that t{uj) grovirs with A'^ as a stretched exponential: 
lnr(c.)~ap2. (7) 

When the initial number of cooperators is not too low, 
the long-lived metastable state is entered prior to fixa- 
tion and the MFT ([7]) is independent of the initial condi- 
tion 12, l^l- Eq. ([7]), confirmed by Fig. [31 implies that for 
scale-free networks with 2 < v < 4, fixation occurs much 
more rapidly than on complete graphs, a phenomenon 
called "hyperfixation" in population genetics [l^ . 

For completeness, we have also studied the class of 
"coordination games" (CGs) characterized (at mean-field 
level) by an unstable interior (coexistence) fixed point 
and two stable absorbing states [H- The fixation proba- 
bility of CGs evolving under the LD on scale-free graphs 
has been found to have the same stretched-exponential 
dependence on N as in the SGs 13|, see Fig. 01 

Discussion & conclusion. We have studied metasta- 
bility and fixation of evolutionary processes on scale-free 
networks in the realm of EGT. Within an individual- 
based approach, we have focused on "snowdrift games" 
with frequency-dependent selection evolving according to 
the LD [10| and characterized by long-lived (metastable) 
coexistence. Exploiting a timescale separation occurring 
at weak selection intensity, we have devised an effec- 
tive (single-coordinate) diffusion theory and, from the 
corresponding backward Fokker-Planck equation, calcu- 
lated the probability and mean time of fixation. These 
quantities exhibit an stretched-exponential dependence 
on the population size, in stark contrast with their non- 
spatial counterparts. Here, important consequences of 
the stretched-exponential behaviors are a drastic reduc- 
tion of the mean fixation time and the possible enhance- 
ment of the fixation probability of a few mutants with re- 
spect to a non-spatial setting. These anomalous fixation 
properties reflect the influence of the network's complex 
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